Bound States in Time-Dependent Quantum Transport: Oscillations and Memory 

Effects in Current and Density 
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The presence of bound states in a nanoscale electronic system attached to two biased, macroscopic 
electrodes is shown to give rise to persistent, non-decaying, localized current oscillations which can 
be much larger than the steady part of the current. The amplitude of these oscillations depends 
on the entire history of the applied potential. The bound-state contribution to the static density is 
history-dependent as well. Moreover, the time-dependent formulation leads to a natural definition 
of the bound-state occupations out of equilibrium. 
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In recent years it has become possible to measure the 
current through single molecules attached to two macro- 
scopic electrodes 1 ' 2 . This is hoped to be a first step to- 
wards the vision of "Molecular Electronics" where single 
molecules become the basic units (transistors, etc.) of 
highly miniaturized electronic devices. 

To address electronic transport on such a small scale 
theoretically, one needs a full quantum description of the 
electronic dynamics. Non-equilibrium Green's functions 
(NEGF) provide a natural framework to study quan- 
tum transport properties of nanoscale devices coupled 
to leads. When a bias is applied, the electrodes re- 
main in local equilibrium while the current is driven by 
the different chemical potentials in the left and right 
lead. In model systems the leads are assumed to be non- 
interacting and the current is computed from the Meir- 
Wingreen formula 3 using approximate many-body self- 
energies Smb- For weakly correlated models £mb ~ 
and the Meir-Wingreen formula reduces to the Landauer- 
Biittiker formula 4,5 , as it should. 

If one wants to account for the full atomistic struc- 
ture of the system, the NEGF formalism is usually 
combined with static density functional theory (DFT) 
and the current is computed from a Landauer-type 
equation&Z&^ i 10 ! 11 ! 12 ! 13 ' 14 . This approach enjoys increas- 
ing popularity, in particular for the description of trans- 
port experiments on single molecules^. From a funda- 
mental point of view, however, the use of static DFT - 
which is an equilibrium theory - is not justified to de- 
scribe non-equilibrium situations. For a critical review 
of this methodology, the reader is referred to Ref. [l5|. 

By construction, the NEGF+DFT approach inherits 
the main assumption of the Landauer formalism that for 
a system driven out of equilibrium by a dc bias, a steady 
current will eventually develop. In other words, the dy- 
namical formation of a steady state does not follow from 
the formalism but rather constitutes an assumption. The 
question of how the system reaches the steady state has 
been addressed theoretically in Refs. Iu^fl7l where it was 
shown that the total current (and the density) reaches 



a steady value if the local density of states is a smooth 
function of energy in the device region. In the same work 
it was also shown that for non-interacting electrons the 
steady current is independent both of the initial condi- 
tions and of the history of the bias, i.e., all memory effects 
are washed out. 

The situation is different, however, if there exist two 
or more localized bound states (BS) in the device re- 
gion, i.e., if the local density of states has sharp peaks 
at certain energies. This case has been investigated in 
Ref. [l8| and further been elaborated in Ref. [l9| where 
it was shown analytically that a non-interacting system 
with bound states exposed to a dc bias does not evolve 
to a steady state. Instead, a steady component of the 
current is superimposed by undamped harmonic current 
oscillations: Assuming that the one-particle Hamiltonian 
H(t) globally converges to a time-independent Hamilto- 
nian H°° when t — > oo, the total current through a plane 
II perpendicular to the transport geometry has the form 
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where the static contribution Ijj is given by the Lan- 
dauer formula. The dynamical part, I^\t), reads 
(atomic units are used throughout) 
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where the summation is over all BS of the final Hamil- 
tonian H°° and the oscillation frequencies are given by 
the BS eigenenergy differences. The quantities A^ b , and 
fb b 1 are defined according to 
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In Eq. Q the double surface integral is over the plane n, 
£ is the embedding self-energy and ^^°(r) are BS eigen- 
functions. The operator /(if ) in Eq. ((4| is the Fermi 
function calculated at the equilibrium Hamiltonian H 
while the state \ip' b ) is related to IV'b ) by a unitary trans- 
formation: \ip' b ) = Mlipif). The "memory operator" M 
depends on the history of the TD perturbation. There- 
fore the coefficients /& &/ are matrix elements of f(H ) 
between history-dependent localized functions. 

In the present work we study these time-dependent 
(TD) phenomena for a number of numerical examples 
using a recently proposed algorithm for TD transport 20 . 
We show numerically that the amplitude of the current 
oscillations may be very large compared to the steady 
component of the current. In addition, the memory 
dependence of the current oscillations will be explicitly 
demonstrated. Furthermore, we address an important as- 
pect which, in the NEGF+DFT approach, has remained 
an open problem, namely how to take BS into account in 
the calculation of the densit y 11 ' 21 . Recently, a somewhat 
empirical scheme was suggested on how to occupy BS 
with energies in the bias window^ 2 -. Here we will show 
how the BS occupations naturally result from the time 
evolution of the underlying KS system. Moreover, for 
non-interacting electrons we provide numerical evidence 
that these occupations (and therefore the density in the 
device region) show a history dependence as well. 

In the long-time limit, the dynamical contribution to 
the density is given by 

n^{r, t) = Y J fb,V cos[(6 fc °° - eg?)t]^°°(r)Vg?(r), (5) 

6,6' 

where the amplitudes fb,v are again given by Eq. ((4]). We 
observe that while the diagonal term, 6 = 6', does not 
contribute to the current of Eq. @ , it does contribute to 
the density of Eq. ([5]) with a history-dependent coefficient 
fb.b- This means that even if we average out the density 
oscillations, history dependent effects will show up in the 
density at the device region. This is a central result of 
our analysis. 

In order to shed more light on the actual dependence 
of the coefficients fb t b' on the history of the TD pertur- 
bation we study model systems of non-interacting elec- 
trons using a recently proposed algorithm 20 suitable to 
describe TD transport through open systems. The cen- 
tral feature of this algorithm is that it allows for the 
numerically exact solution of the TD Schrodinger equa- 
tion in the central device region in the presence of the 
semi-infinite electrodes. 

We consider one-dimensional systems described by the 
TD Hamiltonian 

H(x,t) = -^-^ + U (x) + U(x,t) . (6) 

For times t < the system is in its ground state described 
by the Hamiltonian H°(x) = + U (x). At t = the 
system is driven out of equilibrium by applying a TD 



field U(x,t). We choose the TD perturbation in such a 
way that for t — + oo the Hamiltonian globally converges 
to an asymptotic Hamiltonian, H°° . 

The TD perturbation can be split into three sub- 
regions depending on where in the system it is imposed: 
U a , a = L, R, represents the applied bias in the left (L) 
and right (R) leads and is independent of the position 
within the lead. In the central region, the TD perturba- 
tion (which we call a gate voltage V g (x,tj) may depend 
both on position x and time t. Thus, the total TD per- 
turbation can be written as 

![/l(£) — oo < x < xl 
V g (x,t) x L <x<x R . (7) 
Un(t) Xjt < X < OO 

In the numerical simulations described below, the ex- 
plicit propagation window ranges from xl = —1.2 a.u. 
to xr = 1.2 a.u.. We choose a lattice spacing Ax = 0.012 
a.u. and use a simple three-point discretization for the 
kinetic energy. The initial potential is Uq(x) = for any 
point x in the left or right lead. Therefore, the occupied 
part of the continuous spectrum ranges from momenta 
k = —hp = — v / 2ep to k — kp which is discretized with 
400 fc-points (ef being the Fermi energy). The (non- 
interacting) many-body state is propagated from t = 
to t = 1400 a.u. using a time step At = 0.05 a.u.. 

In the first model the initial state is a Slater determi- 
nant of plane waves with energies less than ef = 0.1 a.u.. 
At t — we suddenly switch on a bias Ul = 0.15 a.u. in 
the left lead and as a result a current flows which after 
some transient time reaches a steady value of about 0.027 
a.u.. After the steady state is attained, at time T = 100 
a.u., we switch on a gate potential 

V g (x,t) = -^^-v g (8) 

for T < t < T + T g with a switching time T g = 20 a.u. 
and the final depth of the potential well v g = 1.3 a.u.. For 
times t > T + T g this gate potential remains unchanged 
at V g (x,t) — —v g and supports two BS at ef = —0.933 
a.u. and = —0.063 a.u.. The resulting TD current 
in the center of the device region is shown in the up- 
per panel of Fig. [1] The development of a steady-state 
current for T < 100 a.u. can clearly be recognized. Af- 
ter the BS are created (t > T + T g ) the current starts 
to oscillate as expected. The amplitude of the current 
oscillation is of the order of 0.35 a.u., i.e., more than 
an order magnitude larger than the steady-state current. 
By Fourier transforming the TD current one can iden- 
tify various transitions which contribute to the oscillat- 
ing behavior. In the long-time limit, only the transition 
between BS survives. In the transient regime, however, 
one finds additional, power-law decaying (1/t) transitions 
from BS to the right continuum at £f as well as to the 
left continuum at ep + Ul—- 

The history dependence of the current oscillations can 
be seen by comparing the current in the upper and lower 
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FIG. 1: Time evolution of the current at x = 0. At t — 
a.u., a bias Ul = 0.15 a.u. is suddenly switched on and the 
system evolves to a steady state. Upper panel: at T = 100 
a.u., a gate potential is turned on (v g = 1.3 a.u. and T g = 20 
a.u.) which creates two BS and results in large amplitude 
oscillations of the current. Lower panel: at T = 100 a.u., a 
first gate potential (v g = 0.2 a.u. and T g = 0) is turned on 
which creates a single BS. Waiting for the transients to de- 
cay, a second gate voltage (v g — 1.1 a.u. and T g = 20 a.u.) is 
then applied which leads to the formation of a second BS and 
therefore to persistent current oscillations. Although H°° is 
identical in both cases, the amplitude of the current oscilla- 
tion is significantly smaller in the second case, illustrating its 
dependence on the history of the system. 



panels of Fig. [T] Both currents were computed by start- 
ing from the same initial state and applying the same 
bias at t = 0. In the lower panel we create the same, fi- 
nal potential as in the upper panel, but in two steps. At 
T = 100 a.u. we suddenly switch on a first gate potential 
with depth v g = 0.2 a.u. which creates one BS. Wait- 
ing for the slow decay of the resulting bound-continuum 
transition we then apply an additional gate potential of 
depth v g — 1.1 a.u. (hence the total depth is 1.3 a.u.) 
with a switching time T g = 20 a.u.. Again we recog- 
nize the persistent current oscillations due to the bound- 
bound transitions. Although the amplitude (about 0.07 
a.u.) is still large compared to the steady-state current, 
it is about a factor of four smaller than the amplitude in 
the previous case. 

Memory effects not only appear in the amplitude of the 
current oscillations but also in the BS contribution to the 
density (Eq. ((SJ)) through the history dependence of 
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FIG. 2: Memory effects for the static part of the density 
in the long-time limit in the presence of BS. The densities 
shown here correspond to systems which are identical to the 
one studied in the upper panel of Fig. [1] except that they are 
computed for three different switching times of the gate po- 
tential of Eq. (0: T g = 0.1 a.u. (solid, black), T g = 5.0 a.u. 
(dashed, blue), and T g = 20 a.u. (dash-dotted, red). The 
inset shows the occupation numbers fb,b (Eq. (j4|) of the two 
bound states as function of switching time. The state with 
lower energy eigenvalue (6 = 1) has higher occupation than 
the one with higher energy. For short switching times the oc- 
cupation is significantly smaller than one, while for adiabatic 
(slow) switching both occupation numbers approach one. 



the coefficients fb,w of Eq. In the long-time limit, 
BS lead to oscillations in (which are connected to 

the current oscillations through the continuity equation) 
and also contribute to the steady part of rS D \ This 
contribution is given by the diagonal part (b = b') of the 
double sum in Eq. ([5]) which implies that fb,b may be 
interpreted as occupation numbers of the BS in the long- 
time limit. In this way, the TD description provides a 
natural way to include BS in a transport calculation. In 
the framework of the Landauer formalism this can only 
be achieved in a somewhat artificial waj*22. 

We emphasize that also the diagonal occupations fb y b 
depend on the history through the memory operator M . 
Here we address the importance of this memory depen- 
dence which is rather crucial in the NEGF+DFT ap- 
proach since BS below the bias window are entirely pop- 
ulated, i.e., fb.b — I- We have computed the time aver- 
aged density n av (x) over an oscillation period for a sys- 
tem with two bound states in the final Hamiltonian. The 
system is the one which leads to the current shown in the 
the upper panel of Fig. Q] except that the gate potential 
is turned on with different switching times T g . In Fig. [2] 
we show n av (x) for three different T g and, as one can see, 
^av differs quite substantially. This difference has to be 
attributed to BS since for the contribution of the scat- 
tering states all memory is washed outii. Of course, the 
relative importance of the BS contribution to n av (x) will 
decrease if e-p (and therefore the contribution of the con- 
tinuum states) increases. By subtracting the continuum 
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contribution from n av and fitting the explicit function 
n^ D ' with the numerical curve (where the only fitting 
parameters are the coefficients fb,b for b = 1,2) we have 
been able to calculate the BS occupations fb,b- Remark- 
ably we have found that fb,b exhibits rather large devia- 
tions from unity, particularly for small switching times 
(see inset of Fig. [5]). As one would intuitively expect, the 
occupation of the state with lower energy (b = 1) is larger 
than for the state with higher energy. In the adiabatic 
limit of very slow switching (T g — ► oo), both occupation 
numbers approach unity which is expected since both 
states are energetically below ep. The occupation num- 
bers also offer an intuitive qualitative picture of the size 
of the current and density oscillations: for relatively short 
switching times (T g < 50 a.u.) the occupation numbers 
deviate substantially from one. Therefore the transition 
probability between the bound states is relatively large 
and so are the oscillations in current and density. On the 
other hand, for large switching times both bound states 
are almost "fully" occupied and the probability of a tran- 
sition between them is small, leading to small amplitudes 
in the dynamical density. 

We have noticed that in the limit of slow switching the 
(transient) transitions between bound states and contin- 
uum can have a rather large amplitude shortly after the 
gate potential is applied, even larger than the (persistent) 
transition between the bound states. As a consequence 



of the very slow decay of the bound-continuum transi- 
tions, one has to propagate for a rather long time before 
the dynamical part of the density takes on the form of 
Eq. ©. 

In summary we have demonstrated that 1) the per- 
sistent current oscillations in the presence of BS can be 
much larger than the steady-state current 2) the ampli- 
tude of the current oscillations can have a strong depen- 
dence on the history of the system 3) a similar history 
dependence is found both for the static and dynamic con- 
tribution of the BS to the density and 4) the occupation 
of the BS is well defined in a TD description of quantum 
transport. 

In the calculations of the present work we assume 
the electrons to be non-interacting. However, the cen- 
tral conclusions (bound-state oscillations in current and 
density, memory effects in both quantities) remain valid 
for any effective single-particle theory such as, e.g., TD 
Hartree-Fock or TDDFT. In this case, the existence of BS 
eigensolution of H°° is in general not consistent with the 
assumption of a steady state and the use of the popular 
NEGF+DFT formalism becomes questionable. 
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